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1. INTRODUCTION 


The Multi-Level Adaptive Technique (MLAT) is a general numerical 
strategy for solving continuous problems such as differential and integral 
equations and functional optimization problems. 

In most numerical procedures for solving such problems, the analyst 
first discretizes the problem, choosing approximating algebraic equations 
on a finite dimensional approximation space, and then devises a numerical 
process to (nearly) solve this huge system of discrete equations. 

Usually, no real interplay is allowed between discretization and solution 
processes. This results in enormous waste: the discretization process, 

being unable to predict the proper resolution and the proper order of 
approximation at each location, uses an approximation space which is too 
fine. The algebraic system thus becomes unnecessarily large in size, while 
accuracy usually remains rather low, since local smoothness of the solution 
is not being properly exploited. On the other hand, the solution process 
fails to take advantage of the fact that the algebraic system to be solved 
does not stand by itself, but is actually an approximation to continuous 
equations, and therefore can itself be approximated by other (much simpler) 
algebraic systems. 

The basic idea of adaptive processes is that an efficient discretization 
of a problem depends on the solution itself: A smooth solution can be 
approximated in a coarse approximation space (a space with relatively few 
degrees of freedom, such as those produced by coarse triangulations). The 
coarse approximation can in fact be very good, provided its order Is 
sufficiently high (e.g., provided finite elements are used which contain all 
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polynomlals up to a sufficiently high degree). A highly-oscillating solution, 
by contrast, can be approximated only in an approximation space which is 
fine enough to resolve the oscillations. In general the solution may be 
smooth in one subdomain, oscillating in another, and may have all kinds of 
singularities around some special points or manifolds. In each region 
then, the efficient discretization will be different in nature. Thus, 
finding the efficient discretization becomes an integral part of the 
problem. The problem is therefore solved iteratively, and at certain stages 
the discretization is adapted to the evolving solution. 

The multi-level techniques go one step further by recognizing that 
it is not necessary, at every stage, to adapt the discretization to the 
solution; it is enough, and much more efficient, to adapt the discretization 
to the error in the solution. A smooth error can very efficiently be 
liquidated by a coarse-space approximation. The fine (and computationally 
expensive) approximation spaces are needed only for approximating the 
highly-oscillating part of the solution; they should be used only to 
smooth out osc i 1 1 at i ng errors. Smoothing the error is much less expensive than 
liquidating it, because it can be done local ly . For example, the process 
of relaxation is a local process which very efficiently smooths the 
error, but which, due to its local character, is very slow in liquidating 
smooth errors. 

Thus, the multi-level adaptive technique Is to use not a single 
but a sequence of approximation spaces (levels), wi th geometrical ly decreasing 
mesh-sizes. New levels may be introduced and changed in the process, and 
they constantly interact with each other. The solution procedure Involves 
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relaxat ion sweeps over each level, coarse-level-to-f ine-level interpolations 
of corrections, and f ine-to-coarse transfers of residuals. This procedure 
has two important basic benefits: On the one hand it acts as a very general 
fast solver of the discrete system of equations. On the other hand It 
provides, in a natural way, a flexible and adaptive discretization. The 
total number n of discrete variables can thus be kept low, and the 
solution of the n algebraic equations is obtained in a low number of 
operations. In fact, the number of operations is only 0(n), since a couple of 
relaxation sweeps at the finest level are enough to make the error so smooth 
that the rest of the work can be done on coarser levels. 

MLAT have so far been developed mainly in finite-difference formulations. 

A comprehensive description of this development, together with historical 

notes, can be found in { B3 1 - for a simplified survey of ideas and software, 

see [B4] . Some further developments are reported in [B5] , [B6] , [B71 • 

and [01]. Throughout that development close agreement has been obtained 

between theoretical predictions and numerical experiments. Both show that 

all boundary-value finite-difference problems considered can be solved by 

the multi-grid algorithm in 5 to 10 work units, where a work unit is 

defined as the computational work involved in processing all the difference equations 

of the finest level one time (i.e., the work in one simple relaxation sweep 

over the finest level). Problems ranged from the simplest model problem, 

through singular pertubation problems, to fairly complicated nonlinear 

Navier-Stokes equations and transonic flows, with a variety of domains and 

boundary conditions. The more recent developments include methods for 

time-dependent problems, very efficient multi-grid embedding processes (e.g., 

for bifurcation problems) and methods for local-truncation extrapolations 

(cf. respectively Sections 3-9, 3.8, and 3.10 below). 



The special capability of the finite-difference multi-level 
structure to create non-uniform, flexible discretization patterns based 
on uniform grids is discussed in Section 3 of [B51 . This capability Is 
obtained by observing that the various uniform grids (levels) need not 
all extend over the entire domain. Finer levels may be confined to 
increasingly smaller subdomains, so as to provide higher resolution only 
where desired. Moreover, one may attach to each of these local Ized finer 
grids’ its own local system of coordinates, to fit curved boundaries or 
to approximate directions of interior Interfaces and thin layers. All 
these patches of local grids interact with each other through the multi- 
grid process, which provides fast solution to their difference equations. 

That structure is very flexible. Since it is based on uniform grids, 
it is feasible to employ high-order difference approximations, and to 
change the order whenever desired, with no extra investment in computational 
work. The effective (the finest) mesh-size can be locally adjusted in 
any desired pattern by changing (extending or contracting) the various 
uniform grids, expending negligible amounts of bookkeeping work and storage. 
These flexibilities of a finite difference technique pose a challenge to 
the finite element method, which traditionally had the edge in terms of 
flexibility, especially in treating curved boundaries. 

It has been theoretically shown in [B31 and (B5) that, using this 
flexible structure in an adaptive solution process governed by certain 
criteria, one can get exponential rates of convergence. That is, the 
error (in solving the differential problen) decreases exponentially as a 
function of the computational work. Such rates are obtainable even for 
problems with singular pertubat i ons , algebraic singularities, corners, etc. 
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Parallel to the above Investigations, an interest has been growing 
in applying multi-level techniques in finite element formulations. Such 
applications are mentioned In {B2] and briefly described In various sections 
of [B31. The purpose of the present article Is to give a full description 
of the multi-level adaptive techniques in terms of finite elements for 
general types of problems, and to discuss their various algorithmic and 

theoretical aspects. 

At the same time other workers, starting with Nicolaides [Nil, ( N2 1 
and then Hackbusch [1111 - [H31 , Bank and Dupont [BDJ and Mansfield [Ml] 
have established the mathematical foundations of certain multi-grid finite 
elements algorithms. For a growing class of problems they have rigorously 
proven the basic multi-grid assertion, namely, that the discrete algebraic 
problem with n degrees of freedom (obtained as a finite element approximation 
to a continuous problem) can be solved in only Cn computer operations; or, 
at least, that Cn log operations are enough to reduce the L norm 
of the error by any desired factor e . The constant C is independent of 
n , but its numerical value Is usually not specified. In fact, if numerical 
values of C were calculated from the rigorous proofs, they would turn out to be 
exceedingly large, much larger for example than practical values of n, 
so that for practical purposes the Cn estimates would look worse than some 
Cn 2 estimates. (The only exceptions are the rigorous estimates in Appendix 
C of [B31 and a similar result in {FI]. But they apply to the mode! 
problem only.) 

Such rigorous investigations are of a very different character than 
the present work. The price of rigor is that the results are far from 
realistic. The proofs give meaningful estimates only for extremely large 
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n, and, even then, the work estimates are orders of magnitude too large. 

The estimates are therefore too crude to yield any practical information, 
e.g., they cannot resolve the difference between more efficient and less 
efficient multi-grid processes. (This difference is crucial in practice; 
it may itself cover several orders of magnitude.) For this reason, and 
since the quantity we try to estimate here is actually nothing but the 
computer time (which of course we know anyway, at least aposteriori) , a 
different type of theoretical studies are preferred by the present 
author . 

Discarding rigor, our studies are based on the observation 
that the important multi-level processes are of a local nature, since low-frequency 
corrections are obtained by coarse-level, processes , wh i ch cost very little. 

One can therefore analyze the crucial aspects of multi-grid processes by 
employing a local mode analysis : Far boundaries and low-frequency changes in 

the coefficients of the equations can be ignored, so that the effect of 
multi-level processes on individual Fourier components can easily be 
calculated. General computer routines have been developed to perform 
this analysis automatically for any given problem, yielding precise 
quantitative predictions of the multi-level efficiency. Experiments with 
various types of equations (see [D1] and [PI]) show the work predictions 
to be precise within a few percent. This tool (combined with related 
observations, the most important of which is the "coupled nonsmoothness" 
mentioned in Section 3-2 below) is therefore useful in selecting efficient 
algorithms (see, e.g., Section 6 and Appendix A in [B3]), in understanding 
the numerical results, and in debugging multi-level programs (see Section 
3.7.9 below). It played an important role in developing to full efficiency 
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many multi-grid finite difference codes. 

The multi-grid experience with finite-difference formulations is 
of course very relevant to finite element ones, since usually finite 
element methods give a type of finite-difference equations. It Is, in 
fact, the only experience we have. Namely, all numerical experiments so 
far, even those based on finite element derivations (see {N31 and [PI]), 
were based on uniform grids (e.g., uniform triangulations), so that 
f ini te - di fference structure of the discrete equations was very explicit. 

But there are some special features ' n finite element formulations which 
do not show clearly in general finite difference equations. (The converse 
is also true.) 

The special featuresof the finite element method show themselves most 
naturally in variational problems, of course. We start therefore our 
study of multi-level finite element methods with the general minimization 
problem, where the entire description is given in terms of the "total energy" 
functional E , the minimization of which is our objective. It is not 
assumed that E is quadratic, nor that it has any other special form, so that 
the described method is applicable to general nonlinear problems. The 
basic multi-level processes (relaxation sweeps, f ine-to-coarse transfers 
of residuals and coarse-to-f ine interpolations of corrections) turn out to 
be completely natural; they are directly determined by the objective 
functional E and the approximation spaces. In particular we reproduce, for 
the general nonlinear case, the observation (made first by Nicolaides, see 
INI]) that the natural residual transfer is the adjoint of the correction 
interpolation. Convergence Is guaranteed since E decreases monotonlcal ly 
by each step of relaxation, as well as by the coarse-level corrections. 

A closer examination reveals that the natural processes are not always 
the best ones. For example, the natural finite-element Interpolation from 



a coarse approximation space to a fine one (e.g., the identity Interpolation 
in case the fine space contains the coarse one) Is sometimes too crude (see 
Section 3.1). In some other cases, Interpolations cruder (i.e., of lower 
order) than the natural one could be used without loss of efficiency (see 
Section 3-3). The i-extrapolatlon technique, which can improve very much 
the nodal values of the approximation (Section 3-10), is better understood 
from a finite-difference point of view. In fact, the extrapolation does not 
considerably improve the quality of the approximation when measured In terms 
of the usual finite element norm. All these findings are related of course 
to the known fact that the point-wise errors of the finite elements are 
quite often much larger than the average error. While this fact has little 
effect in usual finite element algorithms, it is important in the multi-level 
processes, which should take advantage of the relations between coarse and 
fine discretizations. Nevertheless, the natural processes are not bad, and 
can safely be used, even though they may sometimes require an order of 
magnitude more computing time. 

A central issue in finite element methods is how important it is to 
ese uniform subdivisions. The nonuniform elements are important on the 
boundaries, but it seems that in the interior there Is usually no need for 
nonuniformity, while uniformity offers substantial gains in computing time 
and storage. See for example the uniform Interior structure In Figure 2 in 
Section 4.2 below. Uniformity becomes even more advantageous when multi-level 
solution methods are used (one reason being that, since the solution of 
the discrete equations Is much faster, it becomes more Important to speed 
up other parts of the algorithm, in particular the assembly of the stiffness 
equations, which is very laborious when nonuniform elements are used). 
Moreover, the multi-level technique, as mentioned above, has Its own mechanism 
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and a very efficient one, for creating nonuniform discretizations, based 
on a collection of uniform grids. We will discuss this Issue In a later 
part of this paper. 

In the later part of the paper some generalizations will be given. 
Multi-level solutions to general Galerkin (weak form) finite element 
discretizations will be described. This will include new types of relaxation 
(e.g . K distributive relaxation, as mentioned in I B51 , [B6], and £ B71 ) and 
a further study of the relation between interpolations and residual weightings. 
More general problems will be included, such as constrained minimization 
problems, degenerate problems, indefinite and non-elliptic problems, and 
eigenvalue problems. Local mode analysis and adaptation techniques will be 
treated in greater detail. 



- 10 - 


2. MINIMIZATION PROBLEM AMD APPROXIMATION SPACES 

2. 1 The Continuous Problem . 

We consider a d-dlmenslonal variational problem of minimizing a 
functional E(u) over an admissible space S of functions u(x) , where 
x - (xy,...,Xj) and u may be a vector of functions ( u (f) ’* * * ’ u (q)) * ^ 

need not be quadratic, but we assume that F and S are such that a unique 
minimum, denoted* U .exists, at least locally. Thus, our problem Is to 
find U such that 

E(U) - min E(u) . (2 j) 

u€S 

For some standard examples the reader is referred to Chapter U below. 


* Generally, we will use capital letters to denote solutions; the corresponding 

lower-case letters will denote approximations. 
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2.2 Approximation Spaces. 

An approximation space $* Is a finite-dimensional space which 

£ 

contains approximations to U . It is convenient to assume that S 
Is a subspace of S (although this Is not absolutely necessary). The 
dimension of S* is n^ , and its basis functions are <Pj(x), 


(j = i .... ,n^) 


S l - {u* | u*(x) = I u! 4>^(x)} , 

j=, J J 


( 2 . 2 ) 


where uj are the nodal values of the trial function u* . Usually, with 

each approximation space S £ , we associate a real parameter h £ representing 

l 

its typical mesh-size , i.e., the support of each basis function *Pj is 
assumed to have linear dimensions comparable to h^ . 

We may like sometimes to regard S £ as a tensor product of spaces 

S £ ,...,S £ , i.e., 

m sl 

£ I s - 1 .. 1 \ 

u = (u [n , u [2] [q^] ' 

Thus, u £ denotes the j-th nodal-value, u ( £ } is the a-th component- funct 

J J, 2, 

of u £ in S, while u^ is the p-th component-function of u In S . 

The functions Uj £ j need not correspond to u^ . In fact, sometimes 

q^>q , for example, when u £ approximates derivatives of u as well as 

point-values of u itself: 


i on 


Bxi i • • • » 


U (u) 


(x) 


( 2 . 3 ) 


In such cases, we denote a , x»(xj , . . . ,x d ) , v=v 1 + ...+v (J and \£*(vj , . , . ,v d ) 
by a | , xj=(x £ 1 X jd } » V j“ V jl + '*- +V jd 3nd -j“ (v jl V Jd ) ’ 


l l 


respectively. A function u fg] I® formed by EujiPj , summing over 
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£ £ £ 
all J with the same ctj and the same Vj . More general 1y, u [ej 

can be formed by summing ujtpj over allj for which Uj represents the 

same linear combination of expressions like (2.3)* The distinction between 

different functions u,!, Is not essential to multi-level processes, and is 

1 ■ ' 1 l Pi 

not really used below. It can be useful, however, In devising the Interpolation 
routines and in theoretical discussions. 

When u| satisfy (2.3), the approximation elements are called nodal 
finite elements ([SF1, p. 101). 

Multi-level solution processes typically employ a sequence of Increasingly 

I 2 M 

finer approximation spaces S , S , ..., S , wi th corresponding mesh-sizes 
hj > h 2 > ... > h M , and usually h f = 2h Jt+] . S M may be the "target space", 
i.e., the space in which the final approximation is sought, with coarser spaces 

S* , S 2 S M_1 serving only as auxiliary spaces. In adaptive processes, 

however, there is of course no target-space and the relation between spaces 
will become more involved. 

The simplest case is that of the nested approximation spaces 


S 1 c S 2 c ... c S M c S. 


(2.M 


£-1 . £ 

This case may arise if each element of S is a union of S -elements. 

In most applications, however, such a requirement, especially near boundaries, 
would pose a severe limitation. Thus, in our general description we do not 
assume (2.£») , but the reader may like to keep this case in mind o» a simple 
example. 
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2.3 Interpolation Between Approximation Spaces and Relative Smoothness . 

In multi- level processing we will need Interpolation operators 1^ 

k i £ k 

which will be linear transformations from S to S such that l^u Is 

"close" to u . In the nested case S e S It Is natural to take l k 

as the Identity operator I (although this Is not always the best choice). 

Generally, the usual polynomial Interpolation procedures, of suitable order, 

£ k k 

can be employed to obtain nodal values of l k u which approximate u as 

k 

well as possible (i.e., to the order of the best approximation to u from 
S* ). The order of Interpolation may vary by circumstances, and will be 
further discussed below (e.g.. Section 3-1). At this point we only Introduce, 
as a general characterization of 1^ , the matrix I defined by 


,£ k ^ . k£ £ 

'k^i = 1 l ij tp j 
j 


(2.5) 


* u 
Let S K be a space coarser than S . The smoothness of a function 


u € S* relative to S K Is measured by the quantity 


1 


inf 

k c c k 
u € S 


,Jt k l 

L u - u 

II ^ if 


(2.5a) 


Thus, the notion of relative smoothness depends on the interpolation operator 

. For example, in nodal finite elements (2.3), if polynomial interpolation 

k ♦ « 

is applied independently for each function u (gj ‘ re l a t' ve smoothness is 
measured in terms of the smoothness of the individual functions u^ . 

A function u* may be smooth in this sense, without being smooth in the usual 

j? 

sense, e.g., without having smooth component- funct ions u^ . Normally, 

£ 

however, our interpolation operator 1^ will be the natural one, and the 

£ 

smoothness (2.5a) will be equivalent to the smoothness of u ( a ) • 



2 . h Derivatives of Functionals. 

Let E ( ) be any functional defined on S £ . We define Its derlv 1 
atlve with respect to the J-th nodal value by 


E*j (u 1 ) - 1 Im | {E(u £ * 6 ip 1 .) - E(u £ ) } . 

J 6-*0 6 J 


( 2 . 6 ) 


Similarly, we may use the notation E £ j » (E*) £ i and also define derivatives 

k 

with respect to nodal values of another space S : 

E k (u*) = lim -j- fE(u e '+ 5 1^ cp^ ) - £( 11 ^)} . (2.7) 

1 6-0 6 k 1 

If u £ =u k 6 S k cr S l , then E k (u £ ) = E k (u k ) . From (2.5) and (2.7), it follows 
£ 

that If Ej are continuous, then 


rk / t\ r * kt rt 1 ti 

E. (u ) - 1 ! f j A u ) • 


( 2 . 8 ) 


Orders in h. Note that e| is not fully equivalent to the first 
variation (Frechet derivative) of E in S. It is proportional to the 
"volume" |tPj | of the basis function tpj . For example, for nodal finite 
elements (2.3), if u* Is a smooth function, then Ej(u*") = 0(h ). 
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2.5 Rltz Discretizations . 

£ 

As our discrete approximation to (2.1) on S , we have the problem of 
finding U* € S* such that 


E( U*') - min E(u*) . 

l ct .l 

u €S 


(2.9) 


This Implies the "stiffness equations" 
Ej(U*) - 0 . (j-1,2,...,n £ ) 


( 2 . 10 ) 


This is a system of equations for the unknowns (the nodal values 

of u 1 ). 

The (discrete, as well as the original) problem Is called linear If the 
system (2.10) Is a linear system, that Is, If E Is quadratic, or equivalently. 
If the stiffness coefficients E.j(u ) are constants independent of u . The 
system can then be rewritten as 


Z E*.U? = e!(0) 

I J* 1 J 


( 2 . 11 ) 


In multi-level processing, use is made of the residuals of an approx- 

£ no 

imatlon u : the j-th residual is the value of -Ej(u ) , which expresses 

j 

by how much u fails to satisfy the discrete equation (2.10). The discrete 

l 

solution U has zero residuals. 

It is important to notice the different scales of the discrete equations 
(2.10) at different levels l . In finite difference methods we write the 
discrete equations as analogs of the differential equations, and as a result, 
the equations have the same scale on all levels. By this we mean that the 
residuals of any smooth function are all 0(1) in terms of the mesh-size h. This 
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is no longer so In (2.10) where the residuals of smooth functions are 

proportional to |<p|| , which Is often O(h^) • 

In the case of nodal finite elements (2.3), equations (2.10) will 

usually take the form of difference equations. We then denote by m^ 

and m f01 the highest order of differences applied in (2.10) to the 
IP] 

functions U/\ and u * , respectively. Hence, the ratio between the 

la; IpJ m 

residuals of smooth and nonsmooth functions u^ is Ofh^ ) , and for 
U J the ratio is Ofh^ ^) . 

Assembling stiffness matrices for multi-level processing can be 
made by the usual procedures. The need for assembly for several approx- 
imation spaces adds only a small amount to the programming effort and 
costs In computer time and space only a fraction more than the finest-level 
assembly, (cf. Sec. 3 • 7 • 8 - ) 
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3- MULTI-LEVEL SOLUTION PROCESSES 

We consider first the usual situation where some "target" approximation 

space S M is given In which the solution of the stiffness equations (2.10, 

A-M) Is sought. In this section we describe processes which provide fast 

solution to this algebraic system, whether linear or nonlinear, by the 

C 1 c 2 -M-1 

iterative usage of some given coarser approximation spaces b ,b s 

The same processes will later be shown to be the main building blocks In 
more developed adaptive procedures, when no target approximation spaces are 


set in advance. 
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3,1 Coarse* Space Approximation . 

£ 

Let S k be an approximation space coarser than S ; namely, h^ Is 
considerably larger than h £ (typically, h k - 2 h^ ) so that n k <<n t . 
Solving the S k stiffness equations, by any method, Is therefore much less 
expensive than solving the S* equations. The simplest and most fami l 1 ar 
way of using S k in the Iterative solution of the S £ equations Is to Inter 

polate an (approximate) solution from S k to S £ , to serve as a first 

, l 

approximation u : 

i/ ■*- | £ u k , (u k approximates U k ) . (3.1) 


How good this first approximation Is depends on the smoothness of the 

solution U Z . In some cases U* Is so smooth that. If the Interpolation 

I* is of order high enough to exploit the smoothness, then the first 
k 

approximation (3-1) will turn out to be good enough and will require no 

c l 

further Improvement. In such cases, however, the approximation space S 

is not really needed: S k already yields the solution to the required 

accuracy. Thus, if the S £ approximation is at all needed, the first 

approximation (3.1) wl H require a considerable improvement. 

Can we compute a correction to u £ again by the inexpensive use of 

the coarse space S k ? Namely, can we somehow approximate the error 

V £ » U £ - u £ by some V k £ S k ? Normally ', the answer is no. If u k in 

(3.1) is a good enough approximation to U k , then V will be a rapidly 

oscillating function that cannot meaningfully be approximated in the coarse 

l 

space S k . Therefore, before we can reuse coarse spaces, the error V 

*An exception Is the case where S k does not fully use the smoothness of the 
solution U . In that case. If l£ in (3-D is of sufficiently |ilgh order, 
then V* will be smooth enough to be approximated by some V € S . This 
situation is related, however to the inappropriate approximation order being 
used, and will therefore not arise in a fully adaptive procedure. 
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should be smoothed out. Smoothing and coarse-space corrections are described 
In the next sections. Here we add some remarks on the first coarse-space 
approximation (3.1). 

The question arises, of what order and what types should the Interpolation 

(3.1) be7 If S k c S l , a natural choice Is the Identity u* « u k . Usually 

k 

this choice Is not the best: u has an error with large rapid oscillations 

(cf., e.g., [SF] p. 168). Such an error, as we will see below. Is the most 
expensive to liquidate by multi-level processes. Hence, a very substantial 
gain can be made by producing the nodal values of u through a higher-order 
polynomial interpolation, which will give much smaller rapid error oscillations. 
For best results, the rapid error oscillations generated by interpolation should 
not be larger than the rapid (i.e., wavelength 0(h^)) oscillations expected 
anyway in the solution U itself. Hence, for best results, the interpolation 
order should be high enough to exploit all the expected smoothness of the 
solution U. The most significant part of this gain will already be obtained 
if the interpolation order Is such that the order of magnitude of the 
residuals which are produced Is not larger than that of the local truncation 
errors (see Section 3.5 below). In other words, the interpolation should be 
exact for every polynomial which minimizes a functional E both In S and 
in S* (see the example In Section A. 2). In particular, Interpolation of 
nodal values of u^ should be of order m ( a ) + P 0*e*, should use polynomials 
of degree m^+p-1) » w ^ere p Is the order of approximation. For analyses 

of interpolation orders see Section A. 2 In I B3 1 . 
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k k 

Another possibility Is to use two coarse-space solutions, u € S and 

u - 5 € S-J say, with mesh-slzes h.>k k , and to define u* by h-extrapolatlon 

. . J k J 

from u J and u k . This Is a reasonable procedure only when S and S 

have uniform grids. Even then, the accuracy obtained Is at best equivalent to 

that obtained by V , a space with mesh-size hj but with a higher order of 

approximation. In principle, it is less expensive to solve the problem 

than that of S k (since h^ > h R ). Thus, in a fully adaptive process, a 

situation where h-extrapolatlon can profitably be used will not arise. 



- 21 - 


3.2 Error Smoothing by Relaxation . 

A basic solution process, used In multi-level algorithms mainly as an 

error smoothing process. Is the relaxation sweep. The simplest example (and 

In a sense also the best) is the following: 

Gauss-Seldel relaxation sweep : This Is a process In which all the nodal 

£ 

values of some given approximate solution u are scanned one by one In some 

it 

prescribed order. Each nodal value Uj , In Its turn, Is replaced by a new value 

Uj , which Is the value for which the energy (2.10) will be as small as 

possible, other nodal values being held fixed. In other words, u^ Is chosen 

so that the corresponding stiffness equation (2.10) is satisfied. 

I 

Actually, if this equation Is non-linear in u. , it Is better not to 
satisfy it completely, but to make only one Newton step toward Its solution. 
Namely, 


“t t rt / & , J \ / / t , | . 

u. = u. - E . (u ,J ) / E. . (u ) , 
j J J JJ 


(3.2) 


l j l 

where u ,J Is the approximate solution Just prior to replacing Uj by Uj 
(That is, if in our prescribed order i is scanned before j , then 


u. J = u. , otherwise u, - u. .) 
ii i i 

Having completed a relaxation sweep, the system (2.10) is not yet solved, 
of course, because its equations are coupled to each other. A well known and 
extensively used method for solving sparse algebraic systems like (2.10) is by 
a long sequence of relaxation sweeps. This is a convergent process, since 
E(u*) is monotonical ly decreasing. But the rate of convergence is very slow. 
Typically, if n^ is the dimension of S , then the number of sweeps 
required for convergences increases as n^ a (see [ Y 1 ] ) . 

A closer examination, however, reveals that the convergence is not slew 
as long as the error - u^ has rapid oscillations (oscillations with wave 


*,J 
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length comparable to h Jl ) . Such error oscillations are typically reduced 
at least by a factor .5 per sweep. Thus, the convergence slows down only 
when V* becomes smooth. 

In other words, relaxation sweeps, inefficient as they are In solving 

problems, are very efficient In smoothing out the error. This property, which 

will be extensively used below. Is very general for Gauss-Seldel relaxation 

of any non-degenerate discretization of a minimization problem. Degeneracy 

occurs when the stiffness system of equations Is decomposable, at least 

locally, Into several decoupled (or weakly coupled) systems. In such cases 

efficient smoothing can still be obtained by more sophisticated Gauss-Sei del 

relaxation schemes, which take the degeneracy into account, e.g., line 

relaxation In suitable directions. 

\ 

Remark: For a system (q^>l), the relaxation may be slowly converging 

in some , but not all, the components. For example, in the "mixed method", 

when both u,\ and some of its derivatives are taken as independent unknowns, 

(a) 

relaxing over the equations related to derivatives (i.e ( , relaxing 
e!(u*')= 0 over the j for which v ■ w! 0 In (2.3)) will converge very fast. 

This of course does not help, since the lowest-order equations converge slowly; 
but it stresses the fact that these lowest-order equations should be the 
primal concern In the steps below (coarse-grid corrections). These equations 
are also the ones which produce the finite difference analog to the differential 
equations. If there are several zero-order unknowns (uj such that Vj»0) 
per mesh cube, the finite difference analog may be produced only by 
combinations of the corresponding equations. Those special combinations will 
then be slow to converge, so it is no help that some other combinations may 
converge rapidly by suitable relaxation sweeps. 
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Another remark: There are other relaxation schemes which are quite 

natural to the minimization problem. The most obvious one Is the 
steepest descent method, In which all the unknowns are changed simultaneously, 
each one in proportion to the decrease in E per its unit change. In terms 
of difference equations this method is known as Jacobi (under) relaxation . 

The Gauss-Sei del relaxation Is usually preferable, since it requires less 
storage and provides better smoothing rates. Jacobi relaxation, however, 
is more suitable for parallel computations. 

A final remark about relaxation: the efficient smoothing process does 

not continue indefinitely. Except for some ideal cases (e.g., equations with 
constant coefficients in the Infinite domain, uniformly triangularlzed and 
consistently relaxed), a certain level of rapid error oscillations Is always 
coupled to the smooth errors. Starting from a completely smooth error function, 
rapid error oscillations are generated by the relaxation sweeps because of 
boundary interaction and variations In the stiffness coefficients. In particular, 

p 

relaxing with highly oscillatory coefficients will produce ! n rapid error 

m / \ 

oscillations of magnitude 0(h t '“*) times the magnitude of the smooth error. 
This level of " coupled nonsmoothness 1 1 will persist as relaxation slows down. 
Further relaxation sweeps will be wasteful. Moreover, if the error is smoother 
than this level (see footnote in Section 3 . 0 , relaxation may even magnify the 
rapid error oscillations Instead of reducing them. Hence, In such cases it Is 
best to avoid relaxation altogether (cf. analysis in Section A. 2 of [B31). 
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3.3 Measuring Dynamic Residuals . 

In some multi-level algorithms we may wish .to measure the current error 
V 4 » u*’ - u* In order to detect either convergence or slow convergence rates. 

V* cannot be measured dl rectly. Instead, we can measure the residuals -Ej(u ) 
Computing all these residuals Is quite expensive, however. It costs roughly as 
much as a relaxation sweep, so that measuring them after every sweep (as 
requl-red by some algorithms) would double our computational work. Therefore, 
Instead of computing the "static" residuals -Ej(u £ ) , one usually calculates 
the "dynamic" residuals -E^u*’ 1 ) . which are less expensive since they are 
computed anyway In the course of each relaxation sweep (cf. (3-2)). 

It Is Important that the norms in measuring the residuals on different 
spaces S Z are comparable. That Is, they must all be discrete approximations 
to the same continuous norm of the energy first variation. Thus, the L ^ norm 
of the dynamic residuals is given by 


| res. 


n. 

{i 

lj = 1 


[ |4>'f 1 Ej(u*’ J ) ] 2 ] 


1/2 


(cf. Section 2.5), while their norm Is 

|| res. || * » max |q>j | |Ej (u ’^) | . 


(3-3) 


(3.4) 


and the weighted Lj QQrm (perhaps the most useful one) is 


ill'll g ■ j ^ 1 - c(x j ) h « d 


(3.5) 


One (or several) of these norms may be calculated along with each relaxation 


sweep. 
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For algorithmic decisions a dynamic-residuals norm is as good as the static 
one, because (i) fast relaxation convergence must exhibit a fast decrease of the 
dynamic norm; and (II) when the convergence Is slow the dynamic norm Is 
equal to approximately twice the static norm (since the equations are approximately 
symmetric). 
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3.1* Coarse-Space Correction 

As soon as the error V* - U* - u* In S* has been smoothed out by 
relaxation, a good approximation to It can be Inexpensively computed In a 
coarser space. This Idea was used for finite difference equations by 
Southwell [SI] and by others (cf. the historical notes in tB3l), 
and It has a central role In multi-level solution processes. 


I 

Let S k be an approximation space coarser than S 

i.yi 


We wish to wrl te 


equations for V k £ S k so that Its interpolant l*V k '€S- will approximate 
/ as well as possible. Since V* is the solution of the problem 


E(u*+V*) - min E(u*+V*) , 


v l es l 


a natural definition of V k is by the requirement 


E(u Z +I^V k ) 


min Efu^+I^v ) » 

k k K 

V €S 


(3-6) 


which yields the equations (cf. Section 2.*0 


E k ^(u l +lfv k ) - 0 , (I » 1 . ,n k ) • 


(3.7) 

.k 


This is a system of n k equations for the n R nodal values of V 

For general nonlinear problems, or in the case of non-nested spaces 
(S* * S k ) , equations (3-7) are more complicated than necessary. They 
require a special assembly procedure, and the scheme may become complicated 
In later stages(when still coarser spaces are used in solving 0-7) \ see 
Section 3.5). To simplify the scheme we first rewrite (3.7) In the form. 


E kl V + l*V k ) - E k V) - -E k V) 


(3.8) 


Since V k Is smooth, the left-hand side of (3-8) can be approximated by 


E k (l k uV ; ) - E k ( l^u* - ) 


(3.9) 


Observe that we cannot similarly approximate the left-hand side of (3-7) 
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1 k £ 

since I^V may be small compared to the rapidly osci I lating part of u , 

and therefore the error of the coarse approximation may be large compared 

with . In (3.8), by contrast, even If Is smalt, the left-hand 

l k 

side Is still an approximately linear operator in I^V . Indeed, for linear 

I k 

problems and nested spaces (S 3 S') , the left-hand side of ( 3 * 8 ) exactly 
coincides with (3.9). 

Using ( 2 . 8 ) for the right-hand side of ( 3 . 8 ) our new equations become 
:k / 1 k A,,ik\ 

J 


rK/iK X..,.Kv rK/.K _ k£ -It i v 

t.( l £ u +V ) - - I I j j Fj ( u ) . 


( 3 . 10 ) 


If the problem is linear ( E quadratic) these equations can be rewritten as 


a linear system for V 


1 E ia v a " ~ T : ) • 

a J J 


( 3 . 11 ) 


Observe that this problem Is of the same form as the usual S stiffness 

equations ((2.10) or (2.11), for l =* k ), except that the original right-hand 

% 

side Is replaced by linear combinations of residuals from the finer space S 
Thus, no special assembly Is needed for these equations. 

More generally, for nonlinear problems, we introduce the notation 

( 3 . 12 ) 


77k , k l . „k 

U = I „ u + V , 
£ 


in terms of which (3.10) Is written as 

r*W7lk\ nk/.k t- , p ? 7 ^ 

E i (u > ■ M', u ) • E 1 1 j E j<" ) • 


(3.13) 


Again these equations are the usual $ stiffness equations (( 2 . 10 ) for l ■ k ) , 
only the right-hand side Is new, so that no special assembly Is needed. 

k 

The mode of working directly with the coarse-space correction V and 

solving ( 3 . 11 ) Is called the Correction Scheme (CS). The mode of operating 

_k 

with the full approximation (basic approximation plus the correction) U and 
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solvlng (3.13) Is called the Ful 1 -Approx I mat Ion Scheme, (FAS). Note that U , 
„ defined by (3. 12), depends on If l£u l -u‘ .then 5“ coincides 

With U k (the solution of (2.10) for k-l ). At convergence, however, 
u l - U l and V* » 0 , hence 

Tj* 4 - | k |/ (at convergence) • (3.1^) 

£ 

Thus U k Is a coarse-space function which Interpolates a finer-space solution 
and therefore, its nodal values have a finer-space accuracy. This situation 
can be exploited extensively (cf. Section 2 in [B5l) , so that the FAS is not 
only more general than CS, but It also offers other advantages, and Is there- 
fore preferable even in many linear problems* 

Let v k be an approximate solution to (3-11), or u* 4 an approximate 

solution to (3.13), obtained by a method to be specified below (Section 3-6) 

In the first (CS) case, the computed correction should simply be interpolated 
to S and used to correct u , namely 


u* - 
U NEW 


l . .1 k 
U 0LD + ‘k V 


(3.15) 


In the second (FAS) case, it Is Important to realize that u* 4 Itself does not 
approximate a smooth function that can profitably be Interpolated to S . It 
is v k - ^ ~ l k u* , approximating the smooth function V* , which we should 

interpolate. That is, In FAS, 


4w-olo +i ^-'Hld> • 


(3.16) 


Note that \ l \ k Is not the Identity operator, hence (3.16) Is not equivalent 
K & 

to simply u* lju k . This latter Interpolation would destroy the higher- 

H tW K 

e t 

frequency content of u 

The correction (3.15) or (3.16) is called the S k correction to u' . The 


smooth part of the error V practically disappears by such a correct! 


on. High- 
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C, 

frequency error components are Introduced by the l R Interpolation, but 
they are easily liquidated by subsequent relaxation sweeps. Before turning 
to algorithmic details, some remarks should be added concerning the 
Interpolation symbols above. 

Three Interpolation processes were used. The correcti on Interpolation 
(I* In (3.15) or (3,16)), the residual weighting (right-hand side of (3.11) 
and a similar term In (3.13)). and the FAS solution averaging (l fc In (3-13) 
and (3.16)). We use the term "averaging" for Interpolations from a finer 
space to a coarser one, where the coarse space values are obtained as 

weighted averages of fine-space values. 

The correction Interpolation In terms of the above notation is given 


by 


/ . £ k\ — k . k£ 
(1, v ) , - I v, l fJ . . 


’k v 'j 


1 


(3.17) 


It Is thus clear from (3.13) that the residual-weighting Is determined by 
(and actually is the adjoint of) the correction interpolation. For linear problems 
this was observed in [Nl]. We will see later that in more general formulations 
the choice of residual weighting is Independent of the correction Interpolation, 
but In the present case (minimization problems), the relation between the two Is 
natural and need not be violated. Sometimes, if relaxation smooths the residual 
function as efficiently as it smooths V* , simpler and less work-consuming 
residuals weighting can be used without degrading the coarse-space correction 


(see Section 4.2 below), but the possible gain Is quite marginal and 
should not be attempted without specific knowledge. 

The correction interpolation Itself can be made In several ways. In the 

(S* => S k ), it can simply be the Identity. Unlike the first coarse- 


nested case 
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to-flne Interpolation (Section 3.0, higher-order Interpolation Is not 
normally* needed here, since at this stage the smooth error components are 
no longer as dominant. The natural Interpolation (transfer of point-wise 
values of S k ) Is of high enough order to get the full benefit of the 
coarse-space correction. In fact, If high-order elements are used, lower- 
order interpolations could sometimes be used, with no loss in efficiency. 

The Interpolation order need only be high enough so as not to generate 
high-frequency residua,* -arger than the ,ow-f requency residuals of the 
interpolated correction. (See n,ore on interpolation orders In Section A.2 
of U3) and In the chapter on Galerkln formulations.) 

The form of the FAS solution-averaging l( is Immaterial as long as 

o . chr-mirl hp taken when u has wild os ci 1 lations 

u is smooth. Care , however , should be taken wnen 

on the scale of the gr i d , I ,e. , when , ma ^ l u M " u Ml * s com P ara ^ 

I x-y |=h 

to || u*|| . To see this we can view the (nonlinear) equations (2.10) as 
quas i -I I near equations whose coefficients depend on the solution. In the case 
of wild oscillations In u l , the values of l*u' , and hence also the values 
of the coefficients in the coarse-grld equation, (3.|J). depend very much 

on the form of l{ . In such a case the fine-grid difference operator may have 

wildly-oscillating coefficients, and the coarse-grid operator needs to 

represent a proper "homogenization" (cf. Babuska (1975) and Spagnolo (1975)). 

For the purposes of multi-grid processes, enough homogenization Is obtained 

If l k is a "full" averaging operator, i.e., any local operator such that 
l 


j | k w*dx » / w V dx for any w l £ S 


( 3 . 18 ) 


*An exception Is the case mentioned In footnote in Section 3.1, In which 
further gain is obtoined by using higher-order interpolation for the next 


ga 

coarse-space correction. 
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3.5 Relative Local Truncation Errors . 

Let l £ denote an Interpolation (projection) from the solution space 
S to the approximation space S £ . The "local truncation errors" of the 
approximation space S are the values 

7* » e £ (i*u) (3 * 19) 

J J 

which are (up to a sign) the residuals of the true solution U . They are 
so named because they serve as a measure for how well the continuous problem 
is locally approximated by the discrete system (2.10). Approximate knowledge 
of the truncation errors Is important for various algorithmic criteria, such 
as, discretization adaptation criteria and natural stopping criteria. When 
the residuals -eJ(u*) are small In magnitude compared with the corresponding 
truncation errors, then U* Is approximated by u* better than by the exact 
solution U . Hence, at that Instance, U* need no further Improvement, and 
the iterative solution of the stiffness equations may be terminated. 

Moreover, If 7* were known, we could improve the discrete equations. 

In fact, replacing the discrete equations (2 . 10 ) by 

EjVj-tJ. (3 - 20) 

we would get the solution U* = l £ U , that is, a solution which coincides, up 
to the interpolation I* , with the true differential solution U . For 
example, If 1 4 is a point-wi se projection, the nodal values of U* would 

coincide with those of U * 

Of course, since U Is unknown, the truncation errors are not known 

either. Consider, however, the situation described above where we had two 

. £ 0 

approximation spaces, coarser space S and finer one S . Since U is 

usually much closer to U than to U k , we can approximate the S local 
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truncatlon errors Tv by the values 


— k£. rkl |k.,£\ 

t. -■ E,(' a U ) 


(3.21) 


These values are called the local truncation errors of $ relative to 


These values too are not really known until U®" is fully calculated, 
We may, however, replace U 1 in (3-21) by its evolving approximation u . 
More .precisely, IF U* - u* Is a smooth function, then 

E^d^U 1 ) -E*^(I^u £ ) « Ej(U £ ) 

- -E^(u*) 


and hence 


where 


—kl ke, 

T ( «T, 


= i k A\ k V) - ^(u*-) 

I \ i I 


E^(l V) - I l*K(u*) 


( 3 . 22 ) 


? P 


-kfc k?. 


At convergence, when u =U , we indeed get t. * t , 

kf 

The approximate relative-truncation-error r.' is exactly the right-hand 
side of the coarse-space (S^) correction equations (3*13), which may therefore 


be rewritten as 


FjfClft - 


(3.22a) 


We may thus view the role of the fine space (S ) as serving to improve, or -to 
correct, the coarse-space equation ((2.10) for t°k ) by adding to It an 
approximation of the local truncation error. While adding the true truncation 

■ k 

error to the coarse-space equations would make their solution U 

k k 

coincide with the true solution (l.e., U * I U ), adding the relative 
truncation error t. makes U coincide with the fine-space solution 
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(i.e., U k -I k U l ) • 

* — 

T ki< Is a satisfactory approximation to x j as soon as 

|| T kt - < < || T* 1 H • (3^3) 

It has been observed In the numerical experiments that (3.23) Is already 

easily satisfied at the stage of the first transition from the S l relaxation 

sweeps back to S k . Heurlstlcal ly , this is explained as follows: u* at that 

stage' was obtained by the Interpolation (3-1) followed by relaxation. The 

oo — ki 

Interpolation leaves a residuals function E (u ) which Is comparable to t 
both In Its low-frequency components (this is trivial) and In Its high- 
frequency components (this is obtained If a suitable order of Interpolation 
is used, cf. Section 3.1). The relaxation sweeps considerably reduce the high- 
frequency residual components. The bulk of the remaining low-frequency (smooth) 
error In - t k£ is then reduced by adding to T k * the last term in {1.2 

Thus, a relation like (3-23) is satisfied by all components, regardless 
of the norm that is used. In fact, at the said stage, 

£ 

|| x k *- - ^*11 will usually already be comparable to || t || ; further reduction 
of ||T k£ -x kZ || is not meaningful. 

A sequence of refinements . Assume we have a sequence of increasingly finer 
approximations u 1 € S 1 , u 2 € S 2 , . . . , € S* . The corrected coarse-grid 

equations (3.22a) would take the form* 


rk/.-kv k 

E, (U ) - tj 

where T k = if k 

k 

should correct the 5 


(3.24a) 

k 

1-1 . For k<£,-1 , however, the correction Tj 

k+1 

equation not by the original S equation, but by 


*Ve drop the bar from U k . Hereinafter, wherever a flner-space approximation 
u k + ' exists, U k will denote the solution of the correction equation (3-24). 
The original equation (2.10) is the special case . 
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Its own corrected form. Hence, generally (cf. (3.22)) 


(, k+t uk+,) - ] ' u k+1 { E j +1 < uk+1 > - T J + ' } * <" kf (3' 2l|b l 


where 


(3.24c) 


At convergence, the solutions of (3-24) on all levels coincide with each other, 
namely, 

k m |k yk . * b . i 1 iJ’V (at convergence) . 

k+1 k+1 k+2 l 

. k 

Thus, each t* represents the local truncation error on S relative to 

the finest level S . 

Orders in h . By usual Taylor expansions it is easy to find that the 
local truncation error satisfies 


i -l 

Y y 

•tj = b (xj , U) h £ ^ + 0(h Jl J ) , 


(3.25a) 


where the coefficient b depends on the local properties of the exact solution 

, £ £ _ _£ £ 

U but not on the mesh-size h^ , and Yj > Yj (usually Yj = Yj+1 or 

Vf - y*+2 ). yf is the local order of the truncation error. For example, in 

ij ij ij 

the simple nodal case ((2.3), with Vj = 0 and q f ** 1 ) one finds that 
Yj * d+p , where d , the dimension, enters as the scale of the stiffness 
equation (cf. Sec. 2.5) and p is the order of approximation (which usually 
means that the elements contain all polynomials of degree less than p , but 
not all those of degree p ). 

When there are more discrete functions (q 4 >l) then Yj becomes more 

1 k 

complicated, and may actually depend on j . If S and S have the same 


k £ , k l 

structure, and if 1 in S corresponds to J In S (i.e., Xj « Xj 

■ v| ), then they share the same constant In (3.25a), and hence 


and 
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~j * (| + o( ')> • <y ■ y i • y j 1) • <3 - 25b) 

Similarly, we can get a relation of the form 

x^ 1 - c\ l (1 + o(l)) , (3.25c) 

where Cj* is a known quantity (independent of U ). For example, in the 
simple nodal case we have 



(3.25d) 



By combining relaxation sweeps that smooth but the error In u with 

coarse-space corrections that liquidate smooth errors, all error components 

can be efficiently reduced. The question remains of how to solve the 

coarse-grid problem (3-13). This is done by a similar process of relaxation 
k 

sweeps (over u ) combined with sti 1 1 -coarse r-space corrections. 

'Thus, in multi-level solution processes, a sequence of approximation 

1 2 M 1 

spaces S ,S ,...,S is used, starting with a very coarse space S , and endi 

M k 

with the target fine space S . The typical mesh-size of the space S 

is h, «2 ^h . Often, the triangulation for is based on vertical and 

k 0 

horizontal grid lines, and the grid lines of S are every other line of 
„k+1 


Multi-grid algorithms work themselves up from the coarsest level S 

to the finest S M . We will denote by £ the current finest level, that 

k 

Is, the largest k for which an approximate solution u has already been 

£ 

computed. For each i ,a first approximation u is obtained by inter- 

p. i 

oolating from u , and then it is improved by relaxation sweeps and coarse- 
space corrections, using equations (3-24) throughout. One type of multi-level 
algorithm flows as follows (see Figure 1). 

A. Solving on the coarsest grid . Set £=1 . Compute an approximate 
solution u 1 to the stiffness equations (3-24) on the coarsest grid 
(!(«£■ 1) , either by relaxation or by some direct method. (The term direct 
method here means a non-iterative solution of linear systems. If 
the stiffness equations are nonlinear, the direct method will Include a few 
Newton Iterations, where the linear system at each iteration Is solved 
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Figure 1. FAS Multi-Level Algorithm 

The notation -is ge.neAa.Liy exp Pained in the text . 31 jt-/ denote & higheA-oAdea 

{&ee Section 3.1) interpolation, white l£_, ' the ” natuAal " interpolation. 

T k denote* || || . In thi* tfoieehant the coaA*e*t-level (k=l) collection 

equation* aAe *oLved by xelaxation. 
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dtrectly.) Whatever the method, u 1 should be easy to obtain, since S 1 
Is very coarse. 

B. Setting a new finest level . If l -M the algorithm Is terminated. 

If not, Increase i by 1 . Introduce, as the first approximation for the 

new finest level, the function 
i , l l - 1 

u * Vi u 

where* the Interpolation is the higher-order one (see Section 3.1). Having 
assembled the stiffness equations (2.10) for this new level, set k =* £ . 
Generally, k will denote the current operat I on level, (thus, for example, 
when the algorithm later switches to coarser-space corrections, we will have 
k<i ) and u^ € will denote the current approximation on that level. 

Also set sufficiently small, will generally be used as the 

tolerance for solving the k-level equations (3.24). For , a realistic 

value is introduced in Step G below, so the current "sufficiently small" 
value is only temporary. 

C. Starting a new operation level . When we start working on any level 
k , we put * +» , for reasons to become clear in Step E below. 

D. Relaxation sweeps . Improve u by one relaxation sweep (see 
Section 3.2). Concurrently with the sweep .compute some norm e^ of the 
dynamic residuals (see Section 3-3). 

E. Testing the convergence and its rate . If convergence at the current 
operation level has been obtained (e^ $ c^) , go to Step I. If not, and if 
the relaxation convergence rate Is stll 1 satisfactory (l.e., If e. $ rie , 
where n is a prescribed factor which will be discussed below), set 

e^-»-e^ and go back to Step D (continue relaxation). If the convergence rate 
is slow (e^ > ne^ ) , go to Step F (thus obtain coarser-space correction to u^ ) 
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F. Transfer to coarser level . Decrease k by 1. Introduce, as 
the first approximation for the new (the coarser) level k, the function 

J* - 't/*’ (326) 

(see discussion of the FAS solution averaging, Section .J-M. Define the 

S k problem by computing t k as In (3- 24b) , where the first term Is easily 

calculated using (3-26). As the tolerance for this new problem, set 

e -fie , where 5 Is a prescribed factor to be discussed below, 
k k+ 1 

G. Finest level stopping parameter . Concurrently with the computation 
of x k , calculate also Its norm || t k || , using the same norm as used for 
the dynamic residuals (see Step D and Section 3-3). If k = t-1 set 

« # -.||x 1 - , H (3 - 27a) 


Usually we want e £ , the stopping tolerance on the currently finest level, 
to be comparable to || t'|| , so we choose a = (b £ /h f _.j)' r (see (3.25b)), 

or, even more precisely, 


a = 


tv h i/ 

1 - (h t /h k ) p 


(3.27b) 


(cf. (3.25b) and (3 • 25d) ) . 

u Coarse level solution . If k - 1 . solve equation (3.2« by 
relaxation or directly (see Step A above) and 90 to Step I. If k>l , 9° 

to Step C. 

I. Employing a converged so lution to correct a finer level- If k-l , 
go to Step B. If k<t, make the correction (cf. (3-16)) 


k+1 k+1 + |k+1 ( k 

J NEW = U 0LD + k (U 


■ k k+ 1 v 
k+l U 0LD 


(3.28) 


k by 1 and go to Step C. 


and then Increase 
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3.7 Comments 

3.7. 1 Fixed algorithms . The Internal checks In the above algorithm 
can In many cases be replaced by pre-ass Igned flows. (See for example the 
fixed algorithm in Figure 3 in Section 4.2.) In particular, instead of 
checking for slow convergence ( e^ > ne^ , in Step E) , one can simply 
switch to the coarser level k -1 after a pre-asslgned number r £ of 
relaxation sweeps on level k . The parameter r £ , like n , depends on 
the smoothing rate (see below). Hence, If that rate Is known, r may 
be fixed in advance. Similarly, instead of checking for convergence 
(e. £ e. ) , the switch to a finer level k +1 may be made after a total of 

K K 

r^. relaxation sweeps has been made on level k since the last "visit" to 

the finer level. (Thus, if vr *r,<(v+l)r , then v switches from k 

c* c 

to the coarser level k-1 are performed before switching to k+1 .) In 
some cases r ^ a r^ depends on k , and in particular it may have special 
val ues for k = £. . 

The main advantage of fixed algorithms is in saving the work of 
computing the dynamic residuals at each relaxation sweep. This is a significant 
saving when the problem is very simple. For example, in the case of the 
Dirichlet problem (Poisson equation) with linear elements on a uniform square 
grid (Sec. 4.2), a relaxation sweep costs 5 operations (only additions) per node 
A sweep that also computes (3-3) costs 8 operations (7 additions and 1 multi- 
plication) per mode. In more complicated problems, however, the saving Is 
marginal, since a sweep Involves many more operations, and computing the norm 
still adds only 3 operations per nodal value. 

In simple problems, where the fixed algorithms are needed, they are as 
efficient as the "accomodatl ve" ones. In fact, the latter behave like fixed. 
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3 . 7 . 2 . Accumulated work units . For comparison purposes. It Is common 
to record the amount of work performed by a multi -grid algorithm In terms of 
work-units . The work-unit Is the computational work of expressing the stiff- 
ness (finite difference) equations on the finest level M . A relaxation 
sweep on level k usually costs 2 d ^ work units, and so does the 
calculation of the residuals In the transfer from level k to level 
k-1 (except when Injection Is used). All other computational work (mostly 
Interpolations) Is quite small and not easy to express In work units, and Is 
therefore usually neglected. We thus get a measure of the expended work 
which is independent of the hardware and software being used. This measure 
is very convenient in comparing various algorithms and in comparing numerical 
results with theoretical predictions. 

We have examined, both experlmenta 1 ly and theoretically, various types 
of problems (mostly in terms of finite differences). Including systems such 
as the incompressible Navier-Stokes equations. In all cases, the amount of 
work required to solve the problem to the level of truncation errors (3-27) 
was between 5 and 10 work-units. In any case, the number of work-units 
required depends only on the properties of the local opeiator, and not on 
the specific data (forcing terms, boundary shape and boundary conditions). 


3.7.3 Storage requirements . At any given time, the algorithm has at 

k M-k 

most one approximation u on each level k . Assuming h^/h^ ps 2 , 

and denoting by the amount of storage needed for level k (which Is 

usually a fixed multiple of , the dimension of S ), It is clear that 

"n, » 2 d ^ k ' M ^n u , so that the total storaqe is less than 
k M 


(1 +2 _d + 2" 2d + . . . )n 


1 


1 - : 


,-2d M ’ 
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Thls I* only a fraction more than rf M , the storage required for the "target 
space" . 

Moreover, a different variant of the multi-level process requires even 

far less storage. To see this, note that the finer level k+1 Is really 

needed only to provide the correction to the coarser- level equation 

(3.24a). This correction depends only on the local behavior of the solution. 

Hence, we never need the entire finer level , nel ther in core nor in external 

k 

storage. A segment of It Is sufficient for computing tj (at all points 
which are a few mesh-sizes inside that segment). 

3.7.4 Optimizations . The effectiveness of any prescribed multi-grid 

algorithm depends only on local properties (since long error components are 

reduced, very Inexpensively, by coarse-level processes), and can therefore be 

calculated, once and for all, for any type of functional E . either by local mode 

analysis or by numerical experiments. In fact, the numerical experiments 

confirm the mode analysis predictions. We can thus predict which of several 

possible alternatives will give better performance. We can therefore optimize 

our algorithm, including the relaxation method (point-wise or line-wise, marching 

directions, relaxation factors, etc.), the order of interpolations, mesh-size 

ratios (h,/ h k+1 ) , and switching parameters ( n and 6 or r c and r f ). 

The relaxation methods and interpolation orders are discussed In Sections 

3.1, 3.2 and 3.4 above. The mesh-size ratio = 2 ' S c ' ose enou 9*' to 

optimum to warrant Its general use (at least in two and three dimensional problems) 

3.7.5 The Switching parameters . The overall efficiency Is not sensitive 

to the precise choice of 6 and n (Steps E and F above). Quite general, good 

values are 6-0.2 and n = maxu(x) , where u Is the smoothing factor of 

x 

relaxation, l.e., the largest amplification factor (per relaxation sweep) of 
high-frequency error modes. The high-frequency modes are those which are not 
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visible In the coarser space; l.e., their projections on the coarser space 
alias with lower modes. Their ampl iflcatlon factors are eas I 1y computed 
by a local Fourier analysis of the relaxation process (assuming constant 
triangulation, constant coefficients and infinite domain), y is a function 
of x since the triangulation and/or the coefficients may vary over the 
domain. In well-designed relaxation schemes y«0.5 , If not smaller. 

Using the switching parameter n - max y means that relaxation is 
discontinued and coarse-space correction is sought as soon as the amplitudes 
of high-frequency residual components are reduced to the level of the lower- 
frequency amplitudes. Indeed, at this level the smoothing process becomes 
less efficient (see the comments at the end of Section 3.l)»but the error 

U* - U k is already sufficiently smooth. 

If coarse space corrections are not efficient enough, y may always 
be increased and 6 decreased a little, e.g., y may be replaced by y 
and 5 by 5/2 . Theoretically optimal values for n and 5 are discussed 

In Appendix A of l B 3 3 - 


3.7.6 The FAS solution weighting has been discussed above (see Section 
3.4). It is important to emphasize that l k+1 in (3-28), 

(3.26) and (3.24b) should all be identical. A common programming 
error is to have them di ffe rent, at least at some special points. In such 
programs, the coarse-space corrections will deteriorate as soon as their 
magnitude becomes comparable to the difference between the different values 

, ,k k+1 
of l k+1 u . 

3.7,7. A case of avoiding relaxation . If u Is known to be so 
smooth that an interpolation of order higher than m^+p in Step B gives 
a much better first approximation, then that interpolation should not be 
followed by relaxation (see footnotes in Sections 3.1 and 3-2). Go from 



Step B to Step F. A better procedure, of course, would be to employ a higher 
approximation order p . If the smoothness Is not known In advance and a higher 
p poses programming difficulties, a better procedure would be to use t-extra- 
polatlons (see below). 

3-7-8 Programming . The stiffness equations have the same form (3-2*0 

on all levels. Hence, with a suitable data structure, only one relaxation 

routine should be written In which the level number (k) Is a parameter. 

k 

The same Is true for all other basic operations, such as assembling, t 

calculations, and interpolations. An example of such a data structure Is 

exhibited in Appendix B of I B 33 , in Section 4 of [B4] and in the programs of 

k 

[M78]. The routine for calculating t can be produced easily as a combination 
of 3 routine*: a residual calculation routine (RESCAL) , a f i ne-to-coarse 

transfer routine (CTF) and a routine (CORSRES) to compute the coarse-space 
additional term (the first term on the right-hand side of (3.24b)). Both 
RESCAL and CORSRES are trivial modi ficat ions of the relaxation routine (RELAX). 
The interpolation routines, Including CTF, can be written once and for all: 
they depend on the data structure and the types of elements, but not on the 
particular functional E . Thus, the programming for each new problem is 
reduced to the programming of a relaxation routine. A considerable expertise 
is needed, however, to construct a fully efficient relaxation routine (see 
Section 3 in [B2], Section 6 in (85). Lectures 5, 6, 7* 8 in [B6? and [ B 7 3 ) - 
3-7-9 Debugging. By printing out every calculated e^ and || || , 

one gets a nice short summary of each run (cf. Appendix B In [B3l). A typical 
behavior should be exhibited, a deviation from which easily detects most bugs. 
Detailed debugging techniques (for finite difference formulations) are listed 
in Lecture 18 in [B6]. Here we should emphasize only one fundamental 
principle: never settle for any convergence rate slower than (or any work-count 
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larger than) the prediction of the Interior local node analysis. Due to 
the Iterative character of the method, programming (as well as conceptual) 
errors often do not manifest themselves In catastrophic results, but rather 
In considerably slower convergence rates. 



3.8 Nonlinear Problems and Continuatio n (Embeddlngl- 

A basic feature of the above algorithm Is that It has basically the same 
efficiency for nonlinear problems as for linear problems. No linearization Is 
required. A difficulty may arise, however, In problems which have more than 
one local minimum. In such cases. If S 1 ' 1 Is too coarse, the first 
approximation (3.25) may lead to solutions u l converging to a wrong local 
. _t 

minimum in S 

•A cormon way to obtain a first approximation u* close enough to the 
desired minimum is by a continuation (or "embedding") process: the problem, 

including its discretization and its approximate solutions, are driven by some 
parameter y from an easily solvable (e.g. linear) problem toward the desired 
problem, in steps 6y small enough to ensure that the solution to the y-problem 
can serve as a good first approximation in solving the (y+6y) -problem (cf. 

Section 8.2 in [ B 3 3 ) - 

Sometimes the intermediate problems are themselves of interest, i.e., 
they all correspond to interesting possible states of the physical system being 
studied. For example, y may be the total load in a structural mechanics 
problem, or the Reynolds number in a flow problem, or a parameter in terms of 
which the boundary conditions, or the shape of the boundary, are expressed, etc. 
It may then be required to solve the i ntermed i ate-y problems as accurately as 
the final-y problem. When the intermediate problems are not of interest, they 
can be solved to a lower accuracy, using coarser grids. In any case, 
the grids for the intermediate problems cannot be too coarse, lest the first 
approximation obtained from them will not lead to the desl red minimum. In other 
words, the intermediate problems should use a sufficiently fine space S , 
either because this is their desired accuracy, or because solution components 
of wavelengths comparable to h^ are needed to separate between attraction 
regions 11 of different minima. 
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Even though the S* high-frequency components are needed in the 
continuation process, in each 6y step they do not usually change much. 

We can therefore employ the "frozen-T" technique described below, and 
perform most of the 6y steps on very coarse spaces, with only a few 
"visits" to S* . Obtaining a fast approximation by such a continuation 
process will normally require less computational work than the work in 
solving (to a better accuracy) the final problem once its first approximation 


i s g i ven . 
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3,9 Evolution Problems, Optimal-Design Problems, and Frozen-T Techniques . 

We often need to solve not just one isolated problem but a sequence 
of simi lar problems depending on some parameter. For example, we may be 
studying the effect of changing some physical parameters on the "performance" 
of a system, where the performance is measured in terms of the solution U 
of the differential problem. We may want to find for what physical parameters 
the performance is optimal. Or, we may need to solve a sequence of problems 
in a continuation process (Section 3-8). Or, as the most familiar case, 
the parameter may be the time t , and the sequence of solutions describe 
the evolution in time of a physical system. The problem in each time-step, 
more generally, will be to solve some "implicit" part of the evolution 
equations. 

Such a sequence of problems can be handled very efficiently by 
multi-level processes. First, we can use the previous solution (or 
extrapolation from several previous solutions) to obtain a good first approx- 
imation for solving the next problem in the sequence. We will then need 
only one multi-level improvement cycle (involving only a couple of relaxation 
sweeps on each level) to get a satisfactory solution for each problem in 
the sequence. 

Moreover, by a full use of the multi-level structure, one can usually 
do much better. One can exploit the very different rates of change of high- 
frequency and low-frequency components. In parabolic evolution problems, for 
example, high-frequency components converge to a steady state much sooner than 
low-frequency ones do. Hence, after the first few steps, the changes in 
the solution are smooth and can be accurately calculated in increasingly 
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coarser approximation spaces* (This* Incidentally, allows the use of 

largo time-steps without employing Implicit equations at all!) Only once 

In a long time should a step be made In the finest approximation space to 

readjust the high-frequency components. 

Thus, generally, we may like to efficiently solve the next problem In 

our sequence by neglecting the changes in the high - frequency components 

(without neglecting the components themselves). The way to freeze the 

k 

frequencies with wavelengths smaller than O(h^) is to freeze t in 
(3- 24a), i.e., to solve the next problem with the loca I -truncat Ion function 
T k that was previously obtained, thereby limiting the current multi-level 
calculations to the k coarsest spaces S ,...,S . This freezes the 
high-frequency part of the solution but retains its influence on the 
coarse-grid equations. 

Denote by k. the level of calculation (k) In the j-th step; 

k ■* 

that is, x ^ is frozen when solving the j-th problem in the sequence. 

The sequence of levels kj can be selected by monitoring the changes in 
the loca 1 -t runcat ion functions T k . One method is essentially as follows. 
Let T k,J be the function T k at the conclusion of the current step j. 
Thus, T k ’-* - T k, j _1 for, and only for, k>kj . At each step j we update 
the quant i t ies 

6 k « || T k ’ J - t k ’ l(k ’ J) II. i(k,j) = max{i| Kj, k,>k+2> , 

(3.29) 

k 

using the same norm as in Section 3-3 above. That Is, 6j measures the 

total change In t k (owing to calculations at level k+1) since the last 

k+2 k+1 

calculation at level k+2 . Since the last visit to S , t has 



-50- 


remalned unchanged. Had we allowed it to change, its changes would roughly 

be 

5 j k • 2 " y 4 j • (3 - 50) 

where y Is defined In (3.25). When (3.30) exceeds a certain 

tolerance, a new calculation at level k+2 is needed, so kj +1 - k+2 is 

chosdn. 

In the case of evolution problems, the contributions of local errors 
like (3.29) accumulate over time. Instead of <5j , we should then use 

Vi - I At II T k>a - , k ' l(k ’j>|| . (3.31) 

J a-l(k.J) ° 

where At ■ t - t is the time interval related to the a-th step. Hence, 
a a a - 1 

2~ Y 6 k is the estimate of the error contr i buted by the freezing of t . 

J J 

Updating tj +1 costs roughly 2 (2+k_M)d work units, and therefore, the 
condition when to activate such an update should be 

2 -T+d(M-k-2) > x t (3-32) 

where A is the marginal rate of exchange between accuracy and work-units. 

That is, A is a preassigned constant which represents the smallest profit 
rate (the smallest added accuracy per work unit) at which we are still 
willing to invest additional work (see Section 8.1 in [ B3 1 ) . 

Let M, • max k be the finest level used up to time t, In an 

I , (X J 

J a <J 

evolution problem. When (3 • 32) Is satisfied for k-Hj-1, It Implies the 
refinement of our system, l.e., the Introduction of a new finest level 
Mj +1 - Mj +1. Thus, (3.32) serves actually as a discretization-adaptation 
test. Taking norms like (3-31) not over the entire domain but In small 
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subdomalns , we can In fact apply this test in subdomains, thus deciding not 
only when to use a finer level but also where to do it. 

A more extensive description of adaptation techniques and their 
discretization patterns is discussed in later chapters. 



I 
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3.10 T-Extrapolatlon 

When and S* are "similar" spaces we can use Taylor expansions 

- fc -1 

and express the true local truncation function t as a known multiple of 
the computed relative truncation function t - Namely, 

if' 1 = cjrf' 1 * (1 + 0 ( 1 )) (3.33) 

(see (3.25c) - (3.25d) above). Replacing t* 1 by its multiple CjTj 1 

is a trivial and pract I cal ly cost-free modification to the above algorithm, 

which can be conveniently appended to Step G. Such a replacement is called 

l - 1 

jocai-truncat ion extrapolation , or briefly, x-extrapoiat ion . It makes U 

t a 

closer to U (instead of. close to U ), and as a result U will (when 
corrected by u at the next Step I) also be closer to U . In fact, the 
new truncation error agrees with the true truncation error up to higher-order 
terms in h , and therefore the solution with T-extrapolation is equi valent 
to a higher-order approximation . This can of course be true only In the 

l 

sense of approximating the nodal values of I U by those of U; in the 

£ £ 
norm of S , || U -U || is as good as any || u -U || can be. 

Experiments indeed show that if U Is smooth enough, then the T-ext rapolated 

l 

solution u is much closer to the true solution U than is the full solution 

l 

U of the difference equation, namely, 

|| u l - l*U || « || U 4 - l*U || (3.3M 

In fact, if l*U is smooth, (3.3*0 holds already after the first multi-level 

cycle to level l (first Step l). See the heuristic explanation In Section 3*5 

above and the numerical example In IB5). If, on the other hand, l*U has no 

0 

smoothness, then the T-extrapolation will not considerably Improve u . But, 
exactly in this case, the one multi-level cycle is enough to reduce || u* — U*' 
well below || I *U - U*' || since, exactly in this case, t* is not considerably 
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smaller than i 1 ’ 1 , (The size of the residuals after Step B Is roughly 
5 4 ” 1 , and the one cycle easily reduces them well below x l .) So the 
point of T-extrapolatlon Is that, after one cycle^brlnglng the to . ta j 
computational cost typical l y to 5 to 7 work units), it produces an appr o^ 
mat Ion u* which Is guaranteed to be no w orse than U* , the full solution 
of the difference equations , with the nice added feature that any available 
smoothness Is automatically exploited to improve u even further. 

Observe that T-extrapolation is made at level fc-l , the finest level 
at which a non-zero approximation to t is available. The correction is 

automatically carried over to coarser levels via (3 - 2^b) . More precisely, 
(3.24) Is replaced by 


(U k ) 

^k 

= r i 1 

(I<k<l) 

(3.35a) 

T . 
1 

E i" k kt , “ k+, > 

k,k + l { gk+l^k+l) . -k + 1 , 
I j J J 

f 


(1<k<fc-2) 

(3.35b) 

'v £- 1 
T . 

1 


- z if: 1 * 1 eV)] 
j ,j j 

(3.35c) 

T . 

0 . 


(3.35d) 


Richardson extrapolation is the classical form of extrapolating in 
terms of the mesh-size h. It uses the approximate solutions themselves 
( u* and u l ‘\ for example), to produce the extrapolated approximation 
(h u 1 ' 1 - h u l )/(h - h .). It should be pointed out that the t-extrapolatlon 

can be used In many cases where the Richardsoi extrapolation cannot. The 

-l 

T-extrapolat ion employs Taylor expansions of the local , truncation error x , 
while the Richardson extrapolation requires such an expansion for the globa l 
discretization error U* - U , which is more difficult to get. 



4. CONCRETE EXAMPLES 

. i General Linear Elliptic Minimization Problem 

Let fl be any open set In . If v * (v| ,. . . .v^) . Is a 

V 1 v d 

multi-index, D V will denote the differential operator D - D 1 Dj , 

where Dj - 3/3Xj , and |v| will denote i ts order |v| - Vj +. . ,+v^ . 

The space of real functions defined on 0 whose square Is integrable will 
be denoted by L 2 (fl) . This is a HI Ibert space with the scalar product 


(u , v) - J u(x) v(x) dx 
Q : 

1 • 

and the norm ||u|| * (u,u) 4 . The Sobolev space ^(n) is the space 6f 

2 

functions In L (n) whose derivatives of order less than or equal to 
o 

m are also In L (fl) , with the scalar product 


(u,v). ■ I (D v u, D V v) , 

|v|<|4 

and norm Hull . * (u,u) * . For functions u,v € , we define - a continuous 

symmetr I c bilinear functional a(u,v) which is called the strain energy. The 
ellipticity condition is that a(u,v) is positive definite; namely, there 
exists an "ellipticity constant" a>0 such that 


a (u ,u) > 



(4.1) 


Our space of admissible functions will be the space S 3 of functions 

U e H*(a) which satisfy the homogeneous essential boundary condition Bu-0, 
where B Is a bounded linear operator. (The Inhomogeneous condition Bu»g 
can also be used, with obvious changes In the formalism below.) The general 
linear d-dlmensional elliptic minimization problem of order* rfi Is to find 


* See footnote on page 55- 
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U such that 

E(U) » min E(u) , 

U £H* 


(4.2) 


where 

E (u) = -a(u,u) - (f,u) . (^.3) 

2 

It is easy to show that the problem (**.2) - (k,3) has a unique 
solution U , which must satisfy 

a (U ,v) = (f,v) for every v € S . i 1 *- 1 *) 


When integrated by parts, this equation usually yields an equation of the 
form 

(LU,v) » (f,v) for every v 6 S . 


so that U satisfies the differential equation LU * f , at least in a 
weak sense. 

cl c M 

Consider now a sequence of approximation spaces 5 or 

the form (2.2) and typical mesh-sizes . As in S , it is easy 

to see In S l that E(u l ) has a unique minimum U , which satisifles 


*0rder here coincides with Its definition In [SF1, where It is denoted by 

ThiTrder of the associated differential operator L (see (*».5)J is 2m. 
In finite difference formulations, including IB2] - (B71 , the order m-2rf» 
of the differential equation is taken as the order of the problem. 


m 
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i(U*,v*) - (f,v ) 


£ £ 

for every v 65 


(*♦. 6 ) 


The stiffness equations (2.10) here take the form 


a(U l ,iPj) - (f *<Pj ) 


(j*^ * • • • « n o ) 


(4.7) 


where ip* , . . . ,<p l are the basis functions of S* . Equation (4.7) is clearly 
n £ 

equivalent to (4.6). This linear case can also be written as the linear 


system of equations 


4, a Jr u y J ’ 


(j“ 1 » • • » a ) 


(4.8a) 


where (cf. (2.11)) 


l 

a . ~ t . 

jy jy 


/ * l \ 
a (4 > y ,«l> 1 


(4.8b) 


fj " • (4.8c) 

Hence, the multi-grid system of equations (3.24) takes the form 

( 1 <k<Jl) (4.9a) 


_ k M k jk 

* a l6 % f 1 ' 




r k ,.k k+l. . r ,k ,k+l <;k+l _ k.l k+1< 

\ a IB k+1 “ 8 + j 'U <f j * a J Y \ ’* 

(1<k<t-1) (4.9b) 




I ’ 


(k*£) 


(4.9c) 


k k 

where f? denotes the corrected right-hand side fj + Tj, and where the 

ranges of the subscripts are 1<i,P<n k , 1<J,Y<n k+1 * The summat,ons range 

k 

of course over small subsets, since most of the coefficients a.g vanish. 
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The Gauss-Seldel relaxation at level k (see (3*2)) Is given here by 


•l & / t ^-ij f ^ \ / a ^ 

“J ■ “j • ( l a Jt U Y f J >/a JJ 


(*. 10 ) 


If r-extrapolatlon (3.35) Is to be used, equation (4.10b) for k-t-1 should 
be replaced by 


ff' ■ '? r ' * «*■« ■* 

. i if]’-* <fj - i a j Y > 


(*.11) 


The equations above are those of the Full Approximation Scheme (FAS). 
In the linear case, we can also use the Correction Scheme (CS) , in which the 
coarser-space functions are correction functions. Equation (4.9b) is then 
replaced by 

= z l^ ktl (f^ 1 - I a£' ujj + ') . dikit-1) (*.12) 

while (*.9a) and (*.9c) remain the same. Instead of the FAS correction (3.20) 
one should of course use the CS correction 


k+1 _ k+1 i k+1 k 

J NEW “OLD + 'k “ 


(*. 13 ) 



ii.2 The P 1 r 1 ch let Problem with Linear Elements 

The standard example of a linear el 1 Iptle minimization problem is 
the problem (4.2) - (4.3) with the strain energy given by the Dlrichlet 
Integral 

a(u,u) - J (It-) 2 + ... + (l^-) 2 dx ...dx (4.19) 

n 3x i d 

and with the homogeneous Dlrichlet boundary condition u(x)*0 for x € 3ft . 

I 

(techn I ca 1 1 y , this boundary condition is introduced by taking S * Hg » closure 
In H of C”(ft) , where C*(ft) are the infinitely differentiable functions 
which vanish outside a compact subset of ft .) In this case the differential 
operator L in (4.5) turns out to be the Laplace operator 



so that the fliln imi z i ng function U satisfies the Poisson equation with Dlrichlet 
boundary conditions: 

-AU =* f in ft , (4.20a) 

U «* 0 on the boundary 3:. . (4.20b) 

We describe the finite element approximations to the problem in 

1 M 

two dimensions (d»2) . Let the approximation spaces S ,...,S be spaces 

of continuous piecewise linear functions based on triangulations as in 

Figure 2. The triangles touching the boundary may have special forms. 

Such special triangles are called irregular triangles. A grid point serving 

as a vertex of one or more Irregular triangles is called an irregular point. 

it~ 1 

All other grid points are called regular. Notice that S 


is usually 




FIGURE 2 Uniform Triangulations 

The triangulations a > le based on uniformly-spaced horizontal and vertical 
grid tines. The grid tines of S l are every other grid tine of S* \ 
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no t contained In S l , since Its irregular triangles are not unions of 
S l triangles. Generally, the nesting condition S 1 " 1 6 S l would pose 

a severe limitation on the boundary elements of S 8, . 

The basis function Is the continuous piecewise linear function 

£ £ 

which vanishes at all grid points, except for (pj(xj)»1 . It Is easy 
to calculate that, if xj Is a regular grid point, then 

If y*j 


l 
a . 

jy 


i *■ JU 
a(4> Y ,<Pj) 


if y 6 Nj 
otherwise > 


(4.21) 


where nJ is the set of neighbors of j, i.e., y 6 N. if the distance 
between xf and x* is . h, . Hence, at a regular point xf we obtain 

j y 1 

the difference equation 


4llf - Z u 8, 
J y€Nj ^ 


ff . 

J 


(4.22) 


^ 2 2 
The "volume" of a regular basis function tPj is h (its base area is 3h 

and its height is 1), hence f l = (f,4>J) Is h 2 times a weighted average of 

f (x) around xj . Thus (4.22) is h 2 times the usual 5-point approximation 

to the Poisson equation (4.20a). 

Smoothing rate . To calculate the smoothing rate of a Gauss-Seidel 

relaxation sweep applied to (4.22), we change the notation slightly: If 

£ 

xj - (ah, ph), where a and 6 are Integers, then we denote Uj by U^g, and 
ff. by f „ , so that (4.22) takes the form 

J 

K.e ■ V..B - U o.6-1 ■ U a,B+1 ■ ««♦■,» ' f «.6 ' (l| - 23) 


jL . ■ 

Let u denote the value of the approximation u. before the relaxation 
a,e J 
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sweep, and u q 0 Its value after the sweep, 
scanned In a lexicographic order, we can write 
In the form 


If the points (ah,Bh) are 
the relaxation equation (3.2) 


4u 


o,B 


U or1,0 U a,0"1 


- u 


a , 0+1 U a+1 f 0 a » B 


(4 . 2 1 *) 


Let us denote the errors before and after the relaxation sweep by 

u «.8 - Vs “" d Vs ■ Vs - Vs • rMpectWel ’'- Sub,ractln9 

(4.24) from (4.23) we get 


^ V a , B " V a-1.B ’ V a,B-1 " V a,B+I V a+1,B 


0 . 


(4.25) 


We now apply the local mode analysis. Ignoring boundaries, we 


can expand the error in a Fourier Series 

i (0^ + 0 2 B) 


v - I v(0) e 

nS |0|i* 


(4.26) 


where Q = (0^ , ) and where we can take |0| = max 1 | » I 1 ^ ’ s ' nce 

a and B are integers. Similarly, 


/ = 51 v(0) e 

° S |9|<7» 


i (0ja + 0 2 B) 


(4.27) 


Substituting (4.26) and (4.27) into (4.25) we get the amplification factor 
for the 0 component 

I. I I i0 1 ,0 2 I 

(4.28) 


p(e) 


'3 
• > 


10. I0 2 

e + e 


■1 

-18, - e 2 

4 - e - e 

v(e) 



those with |0| > ir/2 , hence 


The high-frequency components are 



» '—2 

is the smoothing factor of the Gauss-Seldel relaxation. That is, each 
relaxation sweep reduces all high-frequency error components by at worst 

the factor . 5 • 

I t is possible to ignore boundaries In this analysis, since we 
only state results concerning high-frequency components. The amplification 
factors for Jow frequencies do depend on the boundaries; (*4.28) for low- 
frequencies holds only in some special cases. We can nevertheless deduce 
from (* 4 . 28 ) that for the lowest frequencies, where |o|«0(h), we get 
u(©)*1-0(h^) . It means that we would need 0(h ) sweeps if we wanted to 

use relaxation to reduce smooth errors by some factor (e.g., by .5)* 

Interpolat ions . If j is a regular point, the natural interpolation 

(2.5) given by the linear elements of Figure 2 is 


,1-1, l 

' 1 j 



if 

n 

♦ — , 
X 

/ £-1 , t" 1 \ 

(x M +h t , x |2 ) 

if 

£ 

X j " 

(x i i 1 . X i2 1 i h t ) 

If 

l 

x . ** 
J 

(xf^+h^, x {2 +h t ) 

if 

l 

X . a 

J 

(x Vl" h l’ X ?2 1 " h t> 


otherwise 


(*4 • 


This linear interpolation Is not good enough for the first interpolation 
(3.1) (Step B in the algorithm). For example, if the solution U Is an^ 
cubic polynomial, it is easy to see that the nodal values of all the discrete 
solutions U* coincide with U . Hence, the first-approximation error 
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u * - u*" Is very large compared with the (vanishing) discretization error 
U* - | £ U . Moreover, the error u l - U 1 Is highly oscillatory, and hence 
more expensive to liquidate by multi-level processes. These difficulties 
do not arise if we replace (4.30) by cubic interpolation, I.e., an 
interpolation which reproduces every cubic polynomial. Note also that if 
we measure the error U* - U In L 2 norm, then it does not vanish even 
if U Is a cubic polynomial; || U £ -U|| , - 0 (h 2 ) notwithstanding the 

L 

"superconvergence" of the nodal values. Hence, from the point of view of 

the L 2 error norm, the natural Interpolation ( ^ - 30) is acceptable. 

At later steps in the algorithm (Step I), where corrections are 

interpolated, the natural interpolation is good enough. Its adjoint can 

.. ,k,k+1 

therefore be used as the res idua 1 -weight i ng in Step F U-e., as l.j 
in (3.24b)). It is interesting to note, however, that a small gain in 
efficiency is obtainable by using "injection", i.e., residual weighting 
as in (3 • 24b) w i th 


|k,k + 1 

U 


k+1 k 
if Xj »x. 


0 otherwise 


(4.31) 


Indeed, Injection is less expensive, because it requires in Step F the 
calculation of only one out of four residuals. At the same time the 
convergence- rate per relaxation work-unit (where the work of Step F is 
not taken into account) happens to be slightly better with injection 
than with the residual weighting (4.30). 

Local mode analysis. In a manner similar to (but more 


involved than) the above analysis of relaxation, one can make a Fourier 
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analysis of the full multilevel cycle, including a fixed number of 
relaxation sweeps and one coarse-grld correction cycle (see 101 ]) . Such 
an analysis shows that with the residual weighting (4.30), and with three 
relaxation sweeps (per cycle) on each level, the multi -grid convergence 
factor (i.e., the factor by which at worst errors are reduced) per relaxation 
work unit is y is .59+. 02. With injection the factor is 

1 - .547 + .015, which happens to be slightly better. The asymptotic convergence 
factors observed In the numerical experiments are precisely In these ranges. 

If we also Include in the work-units count the work of residual 
weighting (see Section 3.7.2) we get the multi-grid convergence factor per 
work unit y . With injection y = y 12/13 « .57, with the residual weighting 


(4.30) y - y 3/4 * -67* 

Incidentally, a good estimate of the convergence factors can be 
obtained from the smoothing factor ii derived above, whose value is always 
easier to calculate. The estimate (explained in Section 6.3 of [B3 1 ) 


.0 _ 

V w V 




(M2) 


which. in the present case yields 


M 


« 0. 5 


0.75 


.595 • 


This estimate of course is independent of the residual weighting (and will 
be approximately realized only when the residual weighting Is good enough) , 
The above convergence factors Imply that In one multi- 
grid cycle the errors are reduced by a factor smaller than 
.14 (.10 in case of injection). This means that one cycle Is actually 
sufficient, since all we have to reduce the error || u* - U* || by is from the 



'< V 
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magnitude of || U 1 ’ 1 - U || (which is its approximate magnitude after Step B) 
to the magnitude of || U* - U || . which is a reduction by a factor not smaller 

than .25 • . 

The one cycle costs less than 

(3 + |)0 + | + ^jV + --- )= T 

work units (cf. Section 3-7.2) in the case of injection (which itself costs 
-J- work-units). Denoting by W H the total number of work units needed to 
solve the S M problem, we get W M = +-j- . Since ^ W M » 

we get = 52/9* 

Thus , the mu i t i - level ' algor i thm solves this problem in 5.8 work units . 

This estimate is confirmed in many numerical experiments. See for 
example (B5) , where results of further improvements by T-extrapola t i on 
are also exhibited. With careful programming, using the correction scheme 
and residual injection, this algorithm requires less than n ^ addi tions 

and no multipl ications , where this count includes all operations (unlike the 
work unit count, which excluded interpolations) and where n^ is the number 
of degrees of freedom in the finest space S (see Section 6.3 in (B3l). 

The flow of the algorithm is shown in Figure 3. The time measurements of 
this algorithm which were made by Craig Poling at (CASE in 1977 are summarized 


in Table 1 . 



' — . Grid"" 

Compute?- — 

X 

33 x 33 

65 x 65 

129 x 129 

Algor 1 thm 

6600 

.028 

cr\ 

OO 

O 

• 

.303 

- 

FMG 


.011 

.048 

.164 

- 

1 Cycle 

■ 

Cyber 175 

.0117 

.0326 

.1085 

- 

FMG 


.0046 

.0161 

.0628 


1 Cycle 

STAR 100 

- 

.0046 

.0109 

.0347 

1 Cycle 


TABLE 1 Time Measurements of the Multi-Level Algo rithms^ 

F MG -tA the Full Mutti-Gfild algorithm mentioned in the. text, and diagrammed -cn 
FiguAe 3. The 1 Cycle algorithm is the same algorithm, but starting ivhen 
the iift U appAoximation is already given in $ M . The cycle include a a 
total 0(5 3 Gauss-Seidel relaxa-tion At veeps on each level, except on the STAR 
100 computeA, wheAe 2 sioeeps Weighted Simultaneous Vis placement (Ace 
[83j) toeAe used instead. Residual weighting ms made by injection. 





1 

* 


« Cubic interpolation 
* Linear interpolation of corrections 


/ 

/ 


* Residual transfer 

= Residual transfer where r- extrapolation 
can be made 


• r relaxation sweeps 


The stage in the algorithm where the error 
is smaller than the discretization error of that level 


Figure 3. Full Multi -Grid Algorithm 
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